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1 Introduction 

Realistic problems in acoustics and aerodynamics are typically time dependent and are often posed 
on unbounded spatial domains. Their numerical solution, of course, requires a bounded domain 
which may be obtained through the introduction of an artificial boundary. To complete the speci- 
fication of the approximate problem, additional boundary conditions must be imposed there. For 
hyperbolic problems in particular, the behavior of the solution over long times can be strongly 
influenced by these additional conditions. 

In the past decade many authors have developed rational approaches to the choice of boundary 
conditions at artificial boundaries. Most have been based on the asymptotic analysis of propagating 
waves. (For a general discussion of asymptotic expansions of waves see Ludwig [13].) There are (at 
least) two natural parameters which can be used to construct the expansion. One is the frequency, 
leading to the theory of geometrical optics. The construction of boundary conditions based on high 
frequency asymptotics was first described by Engquist and Majda [4] and has since been taken up by 
a number of others. Among the advantages of this approach are the following: geometrical optics is 
inherently local, leading always to local boundary conditions; inhomogeneities in the coefficients are 
easily handled; the construction of the expansions is well understood. A disadvantage, however, 
is that the expansion is subject to smooth errors which may not be small. That is, it can be 
difficult to estimate and control the accuracy of these procedures. Engquist and Halpern [3] have 
attempted to rectify this problem by combining the high frequency approximation with one valid 
for low frequencies. It is nonetheless unclear that the correction enhances the accuracy when the 
solution is far from steady-state. 
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A second parameter which may be employed in the asymptotic analysis is the range of prop- 
agation, leading to a far-held expansion. The great advantage of boundary conditions based on 
far-field asymptotics is that their accuracy may be directly controlled either by taking more terms 
in the expansion or by extending the computational domain. A disadvantage is that far-field ex- 
pansions are generally more difficult to compute and require more restrictions on the coefficients of 
the equation. For the wave equation, Friedlander [5] has developed an asymptotic expansion which 
in n dimensions is given by: 


<t> ■ 
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fit* 


■ ■ >9n-i) 

r J+Tr 


(1.0.1) 


Here, r is a radial coordinate and the 0,’s are angular coordinates. In [1], Bayliss and Turkel use 
this expansion to construct accurate boundary conditions. 

In this work we construct expansions analagous to (1.0.1) for more general equations. In 
particular we consider problems which are governed in the far field by: 
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( 1 . 0 . 2 ) 


We assume constant coefficients and that the equation is hyperbolic. We further assume that the 
i-axis is a time-like direction and that all directions in the x plane are space-like. These conditions 
are algebraically equivalent to the following: the (n + l) X (n + l) symmetric matrix C given by: 




1 i = j = 1 
h f 3 - M > 1 

i ’ 

-<*« *,3 > 1 


(1.0.3) 


has one positive and n negative eigenvalues and the n x n matrix A is positive definite. Deviations 
from this form such as inhomogeneities, nonlinearities, forcing terms and scattering bodies may be 
included in the analysis if they are restricted to a compact region in space, as must be the support 
of the initial data. 

The outline of the remainder of the paper is as follows: in section 2 we derive a far-field 
expansion for (1.0.2). In section 3 we use it to derive asymptotic boundary conditions and analyze 
their stability and accuracy. In the final section we present two applications of the theory. The 
first application is to the study of unsteady vortical disturbances to the flow over a flat plate and 
the second is to the solution of the Euler equations for time-dependent flow past a cylinder. More 
detailed versions of aspects of this work will appear in [8] , [9] and [10]. 


2 Asymptotic Expansions 

We now consider the asymptotic analysis of solutions of (1.0.2) at large distances. We find it most 
convenient to work with the Riemann function, S(x,t), which is defined by: 

LS = 0, t > 0, x € R n , 

S = 0, S f = 5(x), 
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t = 0. 


(2.0.4) 

(2.0.5) 


( 2 . 0 . 6 ) 


Here, 6{x) is the Dirac delta function and we assume L is given by: 
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with M < 1. This is the convective wave operator (with subsonic convection). The general 
differential operator appearing in (1.0.2) may be put in this form by a change of spatial coordinates 
only. The restriction to coordinate changes which are independent of time is important for avoiding 
nonuniformities in the resulting expansions. To compute S, we temporarily introduce another 
change of coordinates: 

t = t, (2.0.7) 

yi — xi — Mt, (2.0.8) 

Vi = Xi, i # 1. (2.0.9) 

With this change, L becomes the usual wave operator and the initial conditions for S become: 


5 = 0, Sj-MS Vl =6{y). 


( 2 . 0 . 10 ) 


The first condition, however, implies that S Vl is zero, so that the second reduces to Sj = 6. 
Therefore, S is the well-known Riemann function for the wave equation. In two spatial dimensions 
its expression in cylindrical coordinates is: 


S(r, 

Here, p and s are given by: 
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(1 - M 2 ) 

An asymptotic expansion of S, valid for small, takes the form: 


( 2 . 0 . 12 ) 

(2.0.13) 
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General solutions of (1.0.2) may be expressed in integral form (reference [13]): 
•Z r2x 
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p( z > V’> *l)S[x — y,t — p)2nzdridrpdz. 


(2.0.15) 


Here the domain of integration must include the support of the initial data, the inhomogeneous 
forcing terms, and other perturbations. The integration variables z and tp are polar coordinates in 
y and R and 7 are the polar coordinates of x - y. We introduce: 


T = t 


s(0)r, 


(2.0.16) 
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and note that: 


where e may be expanded: 


t - sR = r - A (0, \}j) - e(x, y), 


OO / \ j 


(2.0.18) 


By a change of integration variables and introduction of (2.0.14) the expression for <f> becomes: 

IJj ^ c ij {0,tP)'^)2nzdzdtl>dv- (2.0.19) 

Jo Jo Jo tt >2 (s(0) + p(0))3>» i+j>i r 


Expanding in e and using (2.0.18) we finally obtain the Friedlander form: 




fj{* - s(0)r,0) 


( 2 . 0 . 20 ) 


A more rigorous exposition of the asymptotic analysis may be found in [8], Some of the conver- 
gence analysis of [13] may also be applicable. We note that the functions /y, j > 1 are determined by 
/o which is called the radiation field. Similar results hold for problems in three spatial dimensions. 
The Riemann function, expressed in spherical coordinates in x space, is: 


S{x,t) = 


2x(r 2 + M 2 t 2 — 2 rMt cos <j> sin 0) 2 


-H(t - (r 2 + M 2 t 2 — 2rMtcos^sin0)?), (2.0.21) 


where H is the Heaviside function. This leads to an expansion for a general solution, u, of the 
form: 

«(*, 0 ~ t ( 2 . 0 . 22 ) 

3-0 r 


— j ((1 ~ -^ 2 (1 ~ sin 2 0 cos 2 <}>))% - Msin0 cos<£). 


3 Asymptotic Boundary Conditions 

Given the expansions (2.0.20), (2.0.22) the derivation and analysis of asymptotic boundary condi- 
tions may essentially be copied from Bayliss and Turkel [1], We introduce the m’th order boundary 
operator, B m : 



„ , nx d (2« - 1). 

B m-I[( d r+s( 9 )dt + 2r 

t= 1 

(3.0.24) 


, 4m+l . 

B m <i> — 0(r 2 ), 

(3.0.25) 

in two dimensions and 

i=l 

(3.0.26) 


B m 4> = 0(r-( 2m+1 >), 

(3.0.27) 


in three dimensions. 
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Standard arguments ([l]) lead to well-posedness and error estimates for a finite domain problem 
with the boundary condition: 

B m v = 0. (3.0.28) 

As we have s(0) > 0 (s(0, 4>) > 0), well-posedness follows directly from a theorem of Higdon [11], 
who in a later work [12] describes techniques for the construction of stable discretizations. As in 
[1], the equation may be used to replace the normal derivatives by tangential derivatives. In the 
numerical experiments which follow only B i is used. 

If we let e(x, t) be the difference between the solution on the unbounded domain and the finite 
domain solution we have: 

Le = 0, t > 0, 2 € H C R n , (3.0.29) 

e = e t = 0, t = 0, (3.0.30) 

B m e = o(f2~( 2m+ ^), x e an, ( 3 . 0 . 31 ) 

where R is a measure of the size of H. The well-posedness implies that for an appropriate Sobolev 
norm: 

IM| m ,n < #m,n||-B m e||o,an. (3.0.32) 

The scaling argument of [1] may be applied to our problem. It yields: 

K m>n = 0(R m+ r). (3.0.33) 

This finally yields the error estimate: 

IMkn= OtiT^+S- 1 )). (3.0.34) 

We note that the curves r = cs, c a constant, are the wavefronts defined by geometrical optics. 
Accordingly, boundary conditions based on high frequency asymptotics are related to those we have 
constructed. The first order boundary condition given by geometrical optics is: 

( «^+~)* = °, (3-0.35) 

for which we have: 

||Se||o, an = 0(i2-( 1+!i ^)). (3.0.36) 

Further corrections, however, will in general differ from those we have computed. 

4 Applications 

4.1 Unsteady Vortical Disturbances Around a Flat Plate 

The first application we consider is to the linearized analysis of unsteady vortical disturbances to 
the flow over a flat plate. The details of this study will appear in Hariharan, Ping and Scott [10]. 
This problem describes the influence of a transverse gust that impinges on an airfoil. The airfoil is 
taken to be a flat plate and the goal is to calculate the unsteady response function. Scott and Atassi 
[15] numerically compute the solution in the frequency domain, where the equations reduce to the 
Helmholtz equation. To calculate the response function, they performed separate calculations for 
each frequency. In contrast [10] deals directly with the time domain where the governing equation 
is a convective wave equation. Therefore the boundary condition proposed here may be directly 
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applied. A typical situation is depicted in Figure 1. The problem is symmetric about the line of 
the plate for uniform flows along the direction of the plate. The governing equations are: 


4>tt + 2Uoo<j> xt + <f> xx = c^vV. 


(4.1.37) 


The boundary conditions along the wake line, on the plate and on the open boundaries respectively 
are 


<S + *-£>* = O' 

d4> 

oy 

*> + jjsjl cos ( # W* + “W* 1 + 2^9) = °- 


(4.1.38) 

(4.1.39) 

(4.1.40) 


Here v is the vertical component of the vortical disturbance and <f> is the velocity potential. In 
addition the response function is defined by: 


• f\J) ’ 


(4.1.41) 


where F(u>) is the Fourier transform of the vortical input, G(oj) is the Fourier transform of the 
unsteady scaled lift given by 

^ = xpcUoo (4.1.42) 

and L is the lift over the plate calculated according to the formula 


j pi,. 


(4.1.43) 


Here T is the boundary of the plate. In this problem the vertical component of the transverse gust 
has the form 

v(x,y ,t) = f(x-t). (4.1.44) 

To compute results that correspond to various frequencies one can choose a Gaussian pulse with 
a wide band width. The results given in Figures 2a, b respectively correspond to solutions obtained 
by solving the above problem numerically. The results are compared to a benchmark cited in [15] 
for Mach numbers .5 and .8 respectively. That there is good agreement is clear. Furthermore, the 
time domain computations are far less costly than those in the frequency domain. For full details 
we refer readers to [10]. 


4.2 Euler Equations 

As a second application we consider the Euler equations of inviscid gas dynamics. Using standard 
notation they are: 

+ V ■ (P&) = °> 


Du - 

p Di = -*» 

De _ „ „ 

p— + p v.u = 0. 


(4.2.45) 

(4.2.46) 

(4.2.47) 
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Figure 2: Comparison between the computed response function of a flat plate airfoil in a transverse 
gust at M = .5 and M = .8 and a Possio Solver (Ref. 15). 
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This system, which must be supplemented by an equation of state relating the density, p, the 
pressure, p, and the internal energy, e, is known to be hyperbolic. 

We assume that in the far field the solution approaches a uniform subsonic flow in the x\ 
direction and consider the linearization of (4.2.45-4.2.47) about the uniform state. We find that 
convenient variables for the analysis are the entropy, S, the vorticity, u>, the divergence of the 
velocity field, q, and the pressure, p. The equations then are: 


dS dS n 

~dt +Uoo dI~ 0, 

(4.2.48) 

duj du 

ar + ““al = 0 ’ 

(4.2.49) 

dq dq 1 _ 2 _ 

T’ + ^oo. 1 V P — 0> 

ut C/ X Poo 

(4.2.50) 

dp dp « _ 

di + Uoo dZ + pooa °° q ~ 0 - 

(4.2.51) 


Clearly the equations for S and a / decouple. Their form implies that these variables should be 
imposed at inflow. A third boundary relation may be obtained from (4.2.50,4.2.51). These obser- 
vations concur with the conclusions of recent papers by Gustafsson [6] and Roe [14]. Combining 
equations (4.2.50) and (4.2.51) yields a convective wave equation for the pressure, to which we 
apply the asymptotic conditions discussed above. 

In summary, we impose, 

(s(0) §}+§;+ ^)iP ~ Poo ) = 0, (4.2.52) 

at each point of the artificial boundary. These are supplemented at inflow by: 

S = Soo, (4.2.53) 

oj = 0. (4.2.54) 

We note that this pressure boundary condition has been previously proposed by Bayliss and Turkel 
[2] for use at a downstream boundary. 

We briefly describe the results of a numerical experiment which demonstrates the accuracy of 
our approach. Details will appear in [9], The (nonlinear) Euler equations were solved in the region 
exterior to a cylinder of radius one with initial conditions modeling an impulsive start. We employed 
a TVD scheme with the superbee flux limiter as part of a two-step Lax-Wendroff procedure (see Yee 
[16]). In Figure 3 the instantaneous density contours are shown for four solutions; two computed 
using the asymptotic boundary conditions and artificial boundary locations of 72 = 5 and 7, one 
computed with the pressure specified and an artificial boundary at 72 = 9 and, for reference, one 
computed on a large domain, 72 = 15. The superiority of the results with the asymptotic condition 
is clear. We suspect that further improvements would be obtained with a nonlinear asymptotic 
analysis as was carried out in the isotropic case in [7]. 
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